import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# --- 1. 参数配置 ---
# 请确保这个路径与你的Excel文件实际存放的路径完全一致
# 注意：Windows路径中的反斜杠'\'最好写成双反斜杠'\\'或正斜杠'/'
file_path = r"D:\Mine in WPS\流致振动\1.xlsx"

# 定义要分析的时间窗口 (开始时间, 结束时间) in seconds
# 我们选择三个有代表性的阶段：爬升期、平稳期、下降期
time_intervals = {
    "Ramp-up": (15, 24),    # 流速接近最大值前的阶段
    "Stable": (25, 34),     # 流速平稳的阶段
    "Ramp-down": (40, 49)    # 流速下降过程中的阶段
}

# --- 2. 加载和预处理数据 ---
try:
    # 读取Excel文件，假设没有表头，第一列为时间，2/3/4列为数据
    df = pd.read_excel(file_path, skiprows=1, names=['Time', 'A', 'B', 'C'])
    # 将 'Time' 列转换为数值类型，以便进行比较
    df['Time'] = pd.to_numeric(df['Time'])
    print("文件加载成功，数据前5行：")
    print(df.head())
except FileNotFoundError:
    print(f"错误：找不到文件 '{file_path}'。请检查路径是否正确。")
    exit()

# 为了更好地显示中文和负号
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


# --- 3. 执行FFT分析并绘图 ---

# 创建一个3x3的子图布局，用于存放9张频谱图
# 每一行代表一个测点 (A, B, C)，每一列代表一个时间阶段
fig, axes = plt.subplots(3, 3, figsize=(18, 12), sharex=True)
fig.suptitle('不同测点在不同阶段的频谱分析 (FFT)', fontsize=20)

# 遍历每个测点
for i, section in enumerate(['A', 'B', 'C']):
    # 遍历每个时间阶段
    for j, (phase, (start_time, end_time)) in enumerate(time_intervals.items()):
        
        # --- 数据切片 ---
        # 筛选出当前时间窗口内的数据
        df_slice = df[(df['Time'] >= start_time) & (df['Time'] <= end_time)]
        
        if df_slice.empty:
            print(f"警告：在 {phase} 阶段 ({start_time}s - {end_time}s) 没有找到数据。")
            continue
        
        # 获取位移数据和时间数据
        displacement = df_slice[section]
        time = df_slice['Time']
        
        # --- FFT计算 ---
        # 移除直流分量（减去均值），使振动信号以0为中心
        y = displacement - displacement.mean()
        
        N = len(y) # 采样点数量
        if N < 2:
            continue
            
        # 计算采样周期 T
        # 注意：这里我们假设采样是均匀的，用平均时间间隔来近似
        sampling_rate = N / (time.iloc[-1] - time.iloc[0])
        T = 1.0 / sampling_rate
        
        # 执行FFT
        yf = np.fft.fft(np.array(y))
        
        # 计算频率轴
        xf = np.fft.fftfreq(N, T)[:N//2]
        
        # 计算振幅谱 (单边谱)
        # 振幅需要乘以 2/N 来获得真实振幅
        amplitude = 2.0/N * np.abs(yf[0:N//2])

        # --- 绘图 ---
        ax = axes[i, j]
        ax.plot(xf, amplitude)
        ax.grid(True)
        ax.set_title(f'截面 {section} - {phase} 阶段')
        ax.set_xlabel('频率 (Hz)')
        ax.set_ylabel('振幅 (mm)')
        ax.set_xlim(0, 15) # 可以调整这个范围，先看0-15Hz的主要频率
        
# 调整子图布局
plt.tight_layout(rect=(0, 0.03, 1, 0.95))
plt.show()

# --- 4. 绘制原始数据图以供参考 ---
plt.figure(figsize=(15, 6))
plt.plot(df['Time'], df['A'], label='截面A', color='black', alpha=0.8)
plt.plot(df['Time'], df['B'], label='截面B', color='red', alpha=0.8)
plt.plot(df['Time'], df['C'], label='截面C', color='blue', alpha=0.8)
plt.title('原始位移数据')
plt.xlabel('时间 (s)')
plt.ylabel('位移 (mm)')
plt.grid(True)
plt.legend()
plt.show()
